Create figures for the clinical paper
All data is coming from tables for the paper
# labels for time points
timeLab = c("Baseline", "Post Run-in", "Week 8")
# set colors
recistCol = c("Non-responder" = "red","Responder" = "blue")
cbrCol = c("CR,PR,SD" = "blue","PD" = "red")
subtypeCol = c("HR+" = "red", "TNBC" = "blue")
timeCol = setNames(c("black", 'grey45', 'grey75'), timeLab)
# set root folder
rootDir = "./"
# Load data patient data
dta_pts = read.csv(paste0(rootDir,"input_tables/table_24patient_info.csv"), row.names = 1)
dta_pts$ID = rownames(dta_pts)
# remove patients with no response
dta_pts = dta_pts %>% filter(Best_Resp != "Unknown")
dim(dta_pts)
## [1] 20 22
table(dta_pts$RECIST)
##
## Non-responder Responder
## 15 5
# load cluster proportions
dat = read.csv(paste0(rootDir,"input_tables/table_per_sample_values.csv"))
props = grep("proportion", colnames(dat), value = T)
#===========================
# add patient info
dat = cbind(dat, dta_pts[dat$patientID,])
Pairwise Wilcoxon test was used for each time point comparing responders vs non-responders. And paired Wilcoxon test was used for baseline vs post run-in and post run-in vs week 8. We took the differences between all time points and then compare those differences between responders vs non-responders to see if there is significant difference in changes after treatment
prop_res = c()
for (k in props)
{
prop_res = rbind(prop_res, c(cluster = gsub("_proportion","",k), runWilcox(k, dat, dta_pts)))
}
knitr::kable(prop_res, "simple",digits = 3)
| cluster | BvsT1 | T1vsT2 | BvsT2 | B_RvsNR | T1_RvsNR | T2_RvsNR | d_T1_B_RvsNR | d_T2_B_RvsNR | d_T2_T1_RvsNR | mean_d_T1_B_R | mean_d_T1_B_NR | mean_d_T2_B_R | mean_d_T2_B_NR | mean_d_T2_T1_R | mean_d_T2_T1_NR |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| B | 0.611 | 0.193 | 0.02 | 0.037 | 0.328 | 0.171 | 0.879 | 0.114 | 0.476 | -1.096 | -0.481 | -7.104 | -2.177 | -6.21 | -1.915 |
| Baso | 0.001 | 1 | 0.432 | 0.792 | 1 | 1 | 0.574 | 0.61 | 0.61 | -0.416 | -0.204 | -0.324 | -0.07 | 0.115 | 0.092 |
| CD4_T | 0.818 | 0.432 | 0.77 | 0.234 | 0.064 | 0.61 | 0.506 | 0.257 | 0.114 | -4.535 | -0.108 | 5.56 | -4.778 | 11.409 | -1.602 |
| CD8_T | 0.378 | 0.232 | 0.275 | 0.721 | 1 | 0.914 | 0.799 | 0.476 | 0.914 | 1.289 | 1.996 | 2.246 | 1.511 | 2.785 | 1.425 |
| cDC1 | 0.185 | 0.275 | 0.275 | 0.598 | 0.752 | 0.476 | 0.205 | 0.762 | 0.352 | 0.058 | -0.455 | 0.056 | -0.002 | -0.029 | 0.243 |
| cDC2 | 0.19 | 0.846 | 0.037 | 0.234 | 0.646 | 1 | 0.574 | 0.257 | 0.762 | -0.461 | -0.101 | -0.802 | -0.386 | -0.283 | 0.007 |
| DNT | 0.587 | 0.906 | 0.636 | 0.562 | 0.959 | 0.61 | 1 | 0.914 | 0.476 | -0.026 | -0.128 | -1.618 | -0.143 | -1.584 | -0.243 |
| DPT | 0.818 | 0.16 | 0.77 | 0.383 | 0.574 | 0.476 | 0.959 | 0.476 | 0.352 | 0.024 | 0.068 | 0.313 | 0.62 | 0.26 | 0.456 |
| GMDSC | 0.698 | 0.695 | 0.695 | 0.37 | 0.646 | 0.61 | 1 | 0.61 | 0.61 | -0.003 | 0.035 | 0 | -0.11 | 0.005 | -0.146 |
| MMDSC_I | 0.404 | 0.232 | 0.625 | 0.959 | 0.879 | 0.762 | 0.721 | 0.476 | 0.61 | 0.126 | 1.346 | 0.319 | -0.27 | 0.139 | -0.941 |
| MMDSC_II | 0.487 | 0.922 | 0.695 | 0.958 | 0.833 | 0.61 | 1 | 0.914 | 0.914 | 0.045 | 0.047 | -0.09 | 0.065 | -0.178 | -0.018 |
| MMDSC_III | 0.159 | 0.432 | 1 | 0.721 | 0.879 | 0.762 | 0.799 | 0.914 | 0.914 | 0.501 | 0.04 | 0.124 | 0.294 | -0.112 | -0.192 |
| Mono_I | 1 | 0.922 | 0.557 | 0.646 | 0.879 | 0.61 | 0.328 | 0.476 | 0.257 | 1.554 | -1.996 | 1.086 | -1.226 | -1.009 | 0.789 |
| Mono_II | 0.842 | 0.906 | 0.064 | 0.833 | 1 | 0.748 | 0.833 | 0.914 | 0.914 | 0.051 | 0.099 | 0.143 | 0.367 | -0.084 | 0.455 |
| Mono_III | 0.459 | 0.193 | 1 | 0.234 | 0.646 | 0.762 | 0.574 | 0.61 | 1 | 0.401 | 7.676 | -3.497 | 1.325 | -6.164 | -3.465 |
| Neu | 0.017 | 0.169 | 1 | 0.195 | 0.799 | 0.914 | 0.027 | 0.067 | 0.61 | 0.047 | -3.725 | 3.14 | -1.235 | 3.081 | -0.1 |
| NK_like | 0.517 | 0.322 | 0.105 | 0.104 | 0.027 | 0.476 | 0.721 | 0.352 | 0.476 | 3.156 | -3.65 | 1.064 | 6.383 | -2.293 | 5.625 |
| pDC | 0.051 | 0.625 | 0.232 | 0.562 | 0.328 | 0.171 | 0.442 | 0.61 | 0.476 | -0.723 | -0.522 | -0.595 | -0.521 | 0.223 | -0.771 |
| UA | 0.89 | 0.232 | 0.275 | 0.154 | 0.328 | 0.762 | 0.959 | 0.171 | 0.038 | 0.006 | 0.064 | -0.024 | 0.354 | -0.069 | 0.3 |
write.csv(prop_res,file = "table_cytof_wilcox_props.csv", row.names = F)
#==========================
# file with marker expression
cytof_dat = read.csv(paste0(rootDir,"input_tables/table_cytof_marker_expression.csv"))
markers = colnames(cytof_dat)[4:ncol(cytof_dat)]
#===========================
# add patient info
cytof_dat = cbind(cytof_dat, dta_pts[cytof_dat$patientID,])
cat("The number of patients with CYTOF data:", length(intersect(cytof_dat$patientID, dta_pts$ID)))
## The number of patients with CYTOF data: 17
# print the number of pairs between time points
mat = printPairs(cytof_dat)
## The number of paired samples between Baseline and Post Run-in : 17
## The number of paired samples between Baseline and Week 8 : 10
## The number of paired samples between Week 8 and Post Run-in : 10
# take only baseline and time 1 values
df_BvT1 <- cytof_dat %>% filter(timePoint %in% c('Baseline','Post Run-in'))
# take only time 1 and time 2 values
df_T1vT2 <- cytof_dat %>% filter(timePoint %in% c('Week 8','Post Run-in'))
Pairwise Wilcoxon test was used for each time point comparing responders vs non-responders. And paired Wilcoxon test was used for baseline vs post run-in and post run-in vs week 8. We took the differences between all time points and then compare those differences between responders vs non-responders to see if there is significant difference in changes after treatment
# Calculate p-values from Wilcoxon test: paired for time points and not paired for response
res = c()
# split by clusters
cytof_dat_byClust = split(cytof_dat, f = factor(cytof_dat$cluster))
for(j in names(cytof_dat_byClust))
{
#
d = cytof_dat_byClust[[j]]
for (k in markers)
{
res = rbind(res, c(marker = k, cluster = j, runWilcox(k, d, dta_pts)))
}
}
res = res[order(res[,'marker']),]
datatable(res)
write.csv(res,file = "table_cytof_wilcox_all_markers.csv", row.names = F)
# make boxplots for markers/clusters that have p-value < 0.05 in d_T1_B_RvsNR
sets = data.frame(res)
sets = sets %>% filter(d_T1_B_RvsNR < 0.05 & cluster%in% c('B','CD4_T','CD8_T','DNT','DPT','NK_like',"cDC1"))
for (i in 1:nrow(sets))
{
k = sets[i,'cluster']
j = sets[i,'marker']
d = cytof_dat %>% filter(cluster == k)
df = cbind(protein = d[,j],
d[,c("patientID", "timePoint", "RECIST","Subtype") ] )
p = ggpaired(df , x="timePoint",y="protein", id = "patientID",
line.size = 0.4, title = paste0(j,". ", k),
line.color = "RECIST",palette = recistCol, xlab="Time") +
facet_wrap("RECIST") +
scale_x_discrete(labels=timeLab)
print(p)
}
# additionally, plot CD103 in cDC1 cluster
k = "cDC1"
j = "CD103"
d = cytof_dat %>% filter(cluster == k)
d = CreateAllFacet(d,"RECIST")
df = cbind(protein = d[,j],
d[,c("patientID", "timePoint", "RECIST","Subtype", "facet") ] )
p = ggpaired(df , x="timePoint",y="protein", id = "patientID",
line.size = 0.4, title = paste0(j,". ", k),
line.color = "RECIST",palette = recistCol, xlab="Time") +
facet_wrap("facet") +
scale_x_discrete(labels=timeLab)
print(p)
#ggpaired(df_BvT1 %>% filter(cluster %in% c("GMDSC")),
for (i in markers)
{
p = ggpaired(df_BvT1, x="timePoint",y=i, id = "patientID", line.color = "RECIST", line.size = 1,
palette = recistCol, title = paste0(i,". Baseline vs Post Run-in"), xlab="Time", ylab="MMI") +
scale_x_discrete(labels=timeLab)+
facet_wrap("cluster", scales = 'free') +
theme(axis.text.x=element_text(angle = 45, hjust = 1))
print(p)
}
for (i in markers)
{
p = ggpaired(df_T1vT2, x="timePoint",y=i, id = "patientID", line.color = "RECIST", line.size = 1,
palette = recistCol, title = paste0(i,". Post Run-in vs Week 8"), xlab="Time", ylab="MMI") +
# scale_x_discrete(labels=timeLab)+
facet_wrap("cluster", scales = 'free') +
theme(axis.text.x=element_text(angle = 45, hjust = 1))
print(p)
}
sessionInfo()
## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] Matrix_1.6-3 scales_1.2.1 gridExtra_2.3 DT_0.30
## [5] data.table_1.14.8 ggrepel_0.9.4 ggpubr_0.6.0 ggplot2_3.4.4
## [9] dplyr_1.1.2
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.7 utf8_1.2.3 generics_0.1.3 tidyr_1.3.0
## [5] rstatix_0.7.2 lattice_0.21-8 digest_0.6.33 magrittr_2.0.3
## [9] evaluate_0.23 grid_4.3.1 fastmap_1.1.1 jsonlite_1.8.7
## [13] backports_1.4.1 purrr_1.0.2 fansi_1.0.4 crosstalk_1.2.0
## [17] jquerylib_0.1.4 abind_1.4-5 cli_3.6.1 rlang_1.1.1
## [21] ellipsis_0.3.2 munsell_0.5.0 withr_2.5.2 cachem_1.0.8
## [25] yaml_2.3.7 tools_4.3.1 ggsignif_0.6.4 colorspace_2.1-0
## [29] broom_1.0.5 vctrs_0.6.3 R6_2.5.1 lifecycle_1.0.4
## [33] car_3.1-2 htmlwidgets_1.6.2 pkgconfig_2.0.3 pillar_1.9.0
## [37] bslib_0.5.1 gtable_0.3.4 glue_1.6.2 Rcpp_1.0.11
## [41] highr_0.10 xfun_0.39 tibble_3.2.1 tidyselect_1.2.0
## [45] rstudioapi_0.15.0 knitr_1.45 farver_2.1.1 htmltools_0.5.5
## [49] labeling_0.4.3 rmarkdown_2.25 carData_3.0-5 compiler_4.3.1